/*  sdfasin.cpp  by K.Tsuru  */
// function ID 3102     DRADIX
/**********************************************************************
SDouble arcsin(x) for |x| <= 1.0
[Algorithm]
A.If x > 1/sqrt(2), using a formura
    arcsin x = pi/2 - arcsin ( sqrt(1-x*x) )
  it can be reduced into |x| < 1/sqrt(2).
B.Pay attention a modified addition theorem
    arcsin(x) = arcsin(y) + arcsin(delta),
  where x, y >=0 and
    delta = x*sqrt(1-y^2) - y*sqrt(1-x^2).
  Because this is an identity, the value of "y" can be chosen arbitrarily
  within the domain.Taking x = y gives delta = 0.

Here we think of the following idea.
1.In the double precision it evaluates x' = asin(x").Where the x" in rhs is
  a rough estimated value of x in double precision if long.
2.In the multi-precision it calculates y = Sin(x').The convergence of series
  Sin(x') is fast.Here an equation x' = arcsin(y) stands up in the multi-precision.
3.Evaluates "delta" in the multi-precision.
4.Using a function AsinSeries(x) which provides the value of arcsin x by series,
     Asin(x) = x' + AsinSeries(delta)
  is calculated.Because the value of "delta" is very small less than 10^(-15),
  we get a fast convergence of series.

The algorithm mentioned above enables us to calculate arcsin x very rapidly.
Most of processing time is used in the evaluation Sin(x').A similar method is
used in the function Log(x).
****************************************************************************/
#ifndef SN_H
#include "sn.h"
#endif

#define usesAsinSeriesInAsin 1  // for debug

/****************************************************************************
An auxiliary function
It receives an approximated value of arcsin x in "approX" and set an exact
value within the present precision in it.
*****************************************************************************/
static void AsinApproX(const SDouble& x, SDouble& approX){
    RealSize C;
    SDouble y, Y;
    uint upPrec = 0;
//If |x| << 1.0,temporally raise up the precition to accurately evaluate the
//value of "x-Y (~= x - approX)" i.e."y" below.
    if(x.NetRdxExp() < 0){
        y = x-approX;
        upPrec = (uint)abs(y.NetRdxExp());
        upPrec = y.ProperUpPrec(upPrec);
        if(upPrec) C.SetEffFig(x.EffFig() + upPrec, C.TEMP_EXTEND);
    }
    Y = Sin(approX);    //Here use much time, about 80% of total one. "SinBS(x)" is called
//  y = x*Sqrt(ONE - Y*Y) + Y*Sqrt(ONE - x*x);  y = (x-Y)*(x+Y)/y;
    y = x*Sqrt(ONE - Y*Y) - Y*Sqrt(ONE - x*x);  //same precision above
    if(upPrec) C.SetEffFig(0);
//Set the fixed value of exponent equal to that of additional value "approX"
//to avoid the evaluation of useless lower figures.
    y = AsinSeries(y, approX.RdxExp());
    approX += y;
}
/*********************************
 r = asin(x) for |x| < sqrt(0.5)
**********************************/
static void AsinForSmallX(const SDouble& x, SDouble& r){
    // If x is too small, UNDERFLOW_ERR in doubleD.
    double dblX = doubleD(x, 0);    //rough estimate in double
#if usesAsinSeriesInAsin
    double ef = double(x.EffFig()), c, absX = fabs(dblX);
    uint xfig = x.Last() - x.First() +1u;

    if(x.NetRdxExp() <= -20) c = -1.0;  // |x| < 1.0e-80
    else if(xfig <= 5u){    // x is short.
        c = 0.75*log10(ef) + log10(absX) - 0.71;
    } else {    // x is long but very small.
        c = ef/64.0 + 0.6*log10(absX) - 0.3;
    }
    // x is short or small, then call AsinSeries()
    if(c < 0){  r = AsinSeries(x);  return; }
#endif
    RealSize C;
//changed on Oct 16, 2000
    uint  efig = x.EffFig(), appFig = (efig < 4000u) ? 20u : 200u;
    uint fig = min(efig, appFig);

    r = asin(dblX);
//An intermediate step is inserted.
    C.SetEffFig(fig);
    AsinApproX(x, r);   //rough estimation in "fig" effective figures
    C.SetEffFig(0);
    if(x.EffFig() > fig) AsinApproX(x, r); //full precision
}
SDouble Asin(const SDouble& x){
    if(x.Sign(3102) == 0) return x;     // x = 0
// check |x| <= 1.0
    SDouble r, pi2;
    if(x.Sign() > 0) r = ONE - x;   // must be r = 1.0 - |x| >= 0
    else             r = ONE + x;
    if(r.Sign() < 0) x.SetError(x.DOMAIN_ERR, "Asin", 3102);    // |x|>1.0

    if(r.Sign() == 0){  // |x| = 1.0
        pi2 = MPi2();                       // pi2 = pi/2
        if(x.Sign() < 0) pi2.ChangeSign();  // pi2 = -pi2
        return pi2;
    }

    double dblX = doubleD(x, 0);

    if(fabs(dblX) >= M_SQRT_2){// |x| >= sqrt(0.5).
    //Using arcsin x = pi/2 - arcsin ( sqrt(1-x*x) )
        pi2 = MPi2();                       // pi2 = pi/2
    //Using the present value of "r", get r = 1 - x*x.
        if(x.Sign() > 0) r *= (ONE + x);
        else             r *= (ONE - x);
        SDouble x1;
        x1 = Sqrt(r);
        AsinForSmallX(x1, r);   //do not use recursive call
        r = pi2 -r;
        if(x.Sign() < 0) r.ChangeSign();    // r = -r;
    } else {
        AsinForSmallX(x, r);
    }
    return r;
}